Viscosity driven instability in rotating relativistic stars 
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We investigate the viscosity driven instability in rotating relativistic stars by means of an iterative 
approach. We focus on polytropic rotating equilibrium stars and impose an m = 2 perturbation 
in the lapse. We vary both the stiffness of the equation of state and the compactness of the star 
to study those effects on the value of the threshold. For a uniformly rotating star, the criterion 
T/W, where T is the rotational kinetic energy and W is the gravitational binding energy, mainly 
depends on the compactness of the star and takes values around 0.13 ~ 0.16, which differ slightly 
from that of Newtonian incompressible stars (~ 0.14). For differentially rotating stars, the critical 
value of T/W is found to span the range 0.17 — 0.25. This is significantly larger than the uniformly 
rotating case with the same compactness of the star. Finally we discuss a possibility of detecting 
gravitational waves from viscosity driven instability with ground-based interferometers. 

PACS numbers: 04.40.Dg, 04.25. Dm, 04.30.Db, 97.10.Kc 



I. INTRODUCTION 

Stars in nature are usually rotating and subject to non- 
axisymmetric rotational instabilities (see , Q , Q or Q 
for recent reviews). An exact treatment of these instabil- 
ities exists only for incompressible equilibrium fluids in 
Newtonian gravity (e.g. [a,IS0)- For these configura- 
tions, global rotational instabilities arise from non-radial 
toroidal modes e lmip (to = ±1, ±2, . . .) when /3 = T/W 
exceeds a certain critical value. Here tp is the azimuthal 
coordinate and T and W are the rotational kinetic and 
gravitational binding energies. In the following we will 
focus on the to = ±2 bar-mode, since it is the fastest 
growing mode when the rotation is sufficiently rapid. 

There exist two different mechanisms and correspond- 
ing timescales for bar-mode instabilities. Uniformly ro- 
tating, incompressible stars in Newtonian theory are sec- 
ularly unstable to bar-mode formation when > /3 sec — 
0.14. This instability can grow in the presence of some 
dissipative mechanism, like viscosity or gravitational ra- 
diation, and the growth time is determined by the dissi- 
pative timescale, which is usually much longer than the 
dynamical timescale of the system. By contrast, a dy- 
namical instability to bar-mode formation sets in when 
P Pdyn — 0.27. This instability is independent of any 
dissipative mechanism, and the growth time is the hy- 
drodynamic timescale of the system. 

There are two representative dissipative mechanisms 
that drive the secular bar mode instability viscosity and 
the gravitational radiation, in the absence of thermal dis- 
sipation. The viscosity driven instability sets in when a 
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mode has a zero-frequency in the frame rotating with the 
star , and the first unstable mode in terms of to is the 
to = 2 bar mode. The quasi-static evolution of the star 
due to viscosity driven instability which varies the cir- 
culation of the star, deforms the Maclaurin spheroid to 
the Jacobi ellipsoid in Newtonian incompressible stars. 
On the other hand, the gravitational radiation induced 
instability HE3 

sets in when the backward going mode 
is dragged forward in the inertial frame (see [il 13. lil Hlj 
for reviews), and the mode in terms of to are all unstable 
when it exceeds a certain to. The quasi-static evolution 
of the star due to gravitational radiation induced insta- 
bility, which varies the angular momentum of the star, 
deforms the Maclaurin spheroid to the Dedekind ellip- 
soid in Newtonian incompressible stars. 

The viscosity driven instability, especially to determine 
the critical value either in an incompressible star or in an 
ellipsoidal equilibrium, has been studied in Newtonian 
gravity , in post Newtonian gravity 0, 0] , and in 
full general relativity 0,0,0 by an ellipsoidal approx- 
imation (e.g. 0) or by an iterative evolution approxima- 
tion (e.g. |15|). and shows that the viscosity drives the 
instability to higher rotation rates /3 sec ^ 0.14 as the con- 
figurations become more compact. There is also a study 
of Newtonian compressible stars to determine the critical 
value of viscosity driven instability [l5f: It is found that 
the star becomes secularly unstable at /3 soc « 0.135±0.02, 
depending on the stiffness of the polytropic equation of 
state. 

The aim of the paper is twofold. One is to investi- 
gate the critical value of viscosity driven instability in 
the compressible stars rotating uniformly. The argument 
that the viscosity driven instability plays a role to de- 
form a star from a Maclaurin spheroid to a Jacobi ellip- 
soid is only true in the absence of internal energy, since 
the total energy has a chance to transfer it to the in- 
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ternal energy without emission from the star (e.g. [l8|L 
Here we assume that the cooling timescale of the star is 
shorter than the thermal heating timescale so that the 
thermal energy generated by viscosity is immediately ra- 
diated away. Therefore the picture of the deformation 
process due to viscosity is quite similar to the case of in- 
compressible stars. The stars are usually considered as 
compressible bodies, and therefore it is worth to take such 
effect into account whether there is a significant change 
on the threshold. In this respect the present work ex- 
tends that of Ref. ^3] to the compressible case. 

The other aim of the present work is to investigate the 
effect of differential rotation on the viscosity driven secu- 
lar bar mode instabilities. For a strong viscosity or strong 
magnetic field circumstances, the star maintains uniform 
rotation. However, in nature the star may rotate differen- 
tially as the Sun. Stellar collapses and mergers may also 
lead to differentially rotating stars (e.g. For the co- 

alescence of binary irrotational neutron stars |20l l2ll |22| . 
the presence of differential rotation may temporarily sta- 
bilize the "hypermassive" remnant, which constructs a 
differential rotation. Therefore it is also worth to take 
differential rotation into account to study viscosity driven 
instabilities in rotating relativistic stars. 

This paper is organized as follows. In Sec.[H]we present 
the basic equations of our treatment of general relativ- 
ity. We discuss our numerical results in Sec. IIIII and 
in Sec. IIVI focusing on the viscosity driven instability 
in uniformly and differentially rotating stars. In Sec. [3 
we briefly summarize our findings. Throughout this pa- 
per, we use the geometrized units with G = c = 1 and 
adopt polar coordinates (r, 9, ip) with the coordinate time 
t. Note that Latin index takes (r, 9, ip). 



II. ITERATIVE EVOLUTION APPROACH TO 
DETERMINE THE THRESHOLD OF VISCOSITY 
DRIVEN INSTABILITY 



Equilibrium configuration of rotating 
relativistic stars 



We briefly introduce our basic approach to construct 
rotating relativistic stars and allow nonaxisymmetric de- 
formation induced by "viscosity" . We introduce a nonax- 
isymmetric line element, when the azimuthal component 
ip is separable in the metric |23| , in spherical coordinates 
(t, r, 9, ip) with a quasi-isotropic gauge (e.g. 0)113) as 

ds 2 = -N 2 dt 2 + A 2 (dr - N r dt) 2 +r 2 A 2 (d9 - N 9 dt) 2 
+r 2 sin 2 9 B 2 (dtp -N^dt) 2 , (2.1) 

where N is the lapse, N r , N e , N v correspond to the shift, 
A and B are the spatial metric functions. In the equilib- 
rium state, we only take the azimuthal component of the 
shift N v , the lapse, and two spatial metric functions into 
account of the metric components. Note that all of them 
are functions of r and 9 only. Once we impose the non- 
axisymmetric perturbation in the lapse, we take all the 



components of shift into account, and relax the depen- 
dence of the functions of (r, 6) to (r, 6, ip) for lapse and 
shift |16|. In our present approach the nonaxisymmetric 
terms have been taken into account at least to the 1/2 
post-Newtonian order. Note that the spatial metric func- 
tions keep the functional dependence of (r, 9) as those 
terms are considered as a higher post-Newtonian order 
which only enters in the order of e a mp x (M/R) ~ 10~ 6 , 
where £ amp is the amplitude of the perturbation we im- 
posed in the lapse, M the gravitational mass, R the cir- 
cumferential radius. It is therefore a good approximation 
to drop the nonaxisymmetric contribution from the spa- 
tial part of the metric. To summarize we compute the 
exact relativistic rotating equilibrium star and perturb 
the geometrical quantities in lapse and in shift, neglect- 
ing the azimuthal perturbation in the spatial metric. 
We also adopt the maximal slicing condition, for which 



the trace of the extrinsic curvature Ka vanishes 



K = j tJ K i:j = 0. 



(2.2) 



The gravitational field equations (Einstein equations) for 
the six unknown functions N, N r , N 9 , N v , A, B are 
written as |16( 

Au = ATrA 2 [E + iP+{E + P)U i U i ]+A 2 K i] K lJ 



-V,j/V> + /3), 
1 



(2.3) 



AN l + -V l Vj-iV 3 = -16irNA 2 {E + P)U l 

+NB- 2 K i: >V J (6{] -v) = J\ (2.4) 
A 2 [rsin6l(iVS-l)] = 16?rr sin 6NA 2 BP, (2.5) 

A 2 C = SnA 2 [P+(E + P)UilP] + '^A 2 K i3 K 13 

-V>VV, (2.6) 

where we introduced the auxiliary functions 

u = ]nN, ( = ln(NA), = lnB. (2.7) 

Note that V, denotes the covariant derivative in terms 
of a flat 3-metric, A = V^V 1 the corresponding Lapla- 
cian, A 2 is the 2-dimensional Laplacian, U l is the spatial 
component of 3- velocity 



A 2 = 



JUL 

dr 2 



Id Id 2 
r dr r 2 89 2 



(2.8) 



N r Af° 1 

U r = , U 9 = , U v = —(SI - N v ).(2.9) 

N N ' N y ' v ' 

In Eqs. (E31, 0, and E and U l are the 

energy density and the 3-velocity, measured by the locally 
non-rotating observer: E = 7 2 [p(l + e) + P)] — P, 7 = 
[l-[A 2 [(U r ) 2 + r 2 (U^) 2 ]+r 2 sin 2 9B 2 (ir?) 2 }}- 1 / 2 , where 
p is the comoving rest-mass density, e the specific internal 
energy, P the pressure. 

The equation for the shift, 



AN* + -VWjNi = J 1 



(2.10) 
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can be further simplified by introducing a vector W» and 
a scalar \ according to 



AWi = Ji, 
The shift can then be computed from 



(2.11) 
(2.12) 



(2.13) 



and will automatically satisfy Eq. H2.10(l . The vector- 
type Poisson equation [Eq. (|2.10() ] for N l has hence been 
reduced to four scalar-type Poisson equations for W % and 
X- 

Let us now discuss the matter part. We treat the mat- 
ter as a perfect fluid, the energy momentum tensor of 
which is 

= p (l + e + — ^ u»u u + Pg» v , (2.14) 

where u M is the fluid 4- velocity. We adopt a T-law equa- 
tion of state in the form 



P = (T - l)pe, 



(2.15) 



where T is the adiabatic index. In the absence of ther- 
mal dissipation, Eq. Ij2.15|l together with the first law of 
thermodynamics imposes a polytropic equation of state 



P = K P 1+1 / n , 



(2.16) 



where n = l/(r — 1) is the polytropic index and n a 
constant. 

The relativistic Euler equation is described in axisym- 
metric stationary spacetime as 



~r~ H" vtuyClj =0 (j = r, I 

h u l 



(2.17) 



where h = (1 + e + P/p) is the specific enthalpy. The 
Bernoulli's equation is derived by integrating the rela- 
tivistic Euler equation [Eq. i|2.17[l ] as (see e.g. |26j l 



H-K + R mt = C, 



(2.18) 



where 



H= =ln/i = ln[l+re], (2.19) 

K= —=]nu t = -u-]nj, (2.20) 
J u l 



and C is a constant of integration. We adopt two types of 
rotation profile, uniform and differential rotation, in this 
paper. For a uniformly rotating star, we simply set the 
rotational energy potential R m t = since ttj = 0. For a 
differentially rotating star, we assume a specific type of 
rotation law as 



ip — A 2 ot (Q c — Q), 



(2.21) 

where Q c is the central angular velocity, A TOt the degree 
of differential rotation, in order to integrate Eq. (|2.17(l 
analytically. In the Newtonian limit (it* — > 1 and u v — » 
w ft, where w is the cylindrical radius), the correspond- 
ing rotational profile reduces to 



A 2 O 

ft = rot c 



H7 2 + A 2 Q 



(2.22) 



which is called ^-constant rotation law. The rotational 
potential energy for this case is 

i? rot = J u'u^dil = ~A 2 Iot (fl c - ft) 2 . (2.23) 

The enthalpy is derived from the Bernoulli's equation 
[Eq. 1(23%!) ] as 

H = H max + (K - K ma3C ) - (Rrot - i^f) 

= ln(l + re max ) + (v max - v) + (ln 7 max - In 7) 
-(R Iot -RZ% (2-24) 

where H max , if max ; i?™ ax represents the values at 
the maximum enthalpy. The rotation law of the star 
[Eq. (J22TJ] becomes 



A ™tWc-n)- (2.25) 



We also rescale the gravitational constant during the 
iteration process in order to determine the radius of the 
star [23| . Since the gravitational constant only enters 
through the matter of the star, we split the equation for 
lapse into two parts: One comes from the spacetime ge- 
ometry and the other from the matter, and solve them 
independently (y — v q + f m , is the contribution from 
the spacetime geometry and u m the one from the matter). 
After that we set the equatorial surface in the computa- 
tional domain and varies the gravitational constant |27^ 
in the following: 



(7J max + ^ ax - ln 7 max + i?™ ax ) - (H sur + z^ ur - ln 7 sur + i?™ r ) 

^ - .siit - .nna.x 1 (2.26) 
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where iT nax , z/™ ax , ln7 max , i? max denotes the value at 
the maximum enthalpy, H sur , z^ ur , In 7 s " 1 ', i? sur denotes 
the value at the equatorial surface of the star, which are 
unknown in each iteration step. 



In order to test our numerical code internally, we check 
the virial identities GRV2 [H and GRV3 (H, the latter 
being a relativistic version of the classical virial theorem. 
The relative errors are defined by 



fn d6 f„°° rdr <7GRV2m (r, 9) 
GRV2 = 1 + Jo J ° GRV2m ; ; , (2.27) 

Jo d0 JO rdr cr GRV2q(f, 0) 

GRV3 = IT dL P Jo sin 0de JiT r2 dr WGRVSqjr, 0, y) + crGRV3m (r, 9, ip)] ^ 
/ * dip sin 0d9 J °° r 2 dr cr Gm3m (r, 6, tp) 



where 



introduce the following nondimensional variables 



a G RV2m(r-,0) = 8irA 2 [P + [E + P)UiU% (2.29) 
<JGW2 q (r,9) = ^A 2 KijK^ -V^VV, (2.30) 



^GRV3m(r,0,^) = 4ttA 2 B[3P + (E + P)UilP 



^GRV3q(r, 0, f) 



-.l 2 A,,/\'< - Vii/VV 



1 



(2.31) 



(2.32) 



Note that the two quantities GRV2 and GRV3, that 
should be identically in the ideal equilibrium config- 
uration, have already been rescaled by the typical source 
term of the equation, and therefore automatically defined 
as a relative error of our computation. 

The gravitational mass M, proper mass M p , total an- 
gular momentum J, rotational kinetic energy T, gravita- 
tional binding energy W can be computed from 

/>2lT />7T />oo 

M= dip sin0<£0 / r 2 dr A 2 B[N[E + 3P 
Jo Jo Jo 

+(E + P)U l U l ] + 2r sin 9B(E + P)N v UiU l ], 

(2.33) 

M p = I dip sin0d9 / r 2 dr A 2 Bjp, (2.34) 
Jo Jo Jo 



r 2 dr 



J = dip sin 

Jo Jo Jo 

xrsin 6A 2 B 2 (E + P)UiU\ 

T = - d(p sin9d9 / r 2 dr 
2 Jo Jo Jo 

xrsm6nA 2 B 2 (E + P)UiU l , 
W = M p + T-M. 



(2.35) 



(2.36) 
(2.37) 



Since we use a polytropic equation of state in the equi- 
librium, it is convenient to rescale all quantities with 
respect to k. Since k"/ 2 has dimensions of length, we 



M = «r"/ 2 M, R = K - n / 2 R, J = n~ n J, 
f = K- n ' 2 T, W = K - n / 2 W, Cl = n n l 2 VL. 



(2.38) 



Henceforth, we adopt nondimensional quantities, but 
omit the bars for convenience (equivalently, we set k = 

Our computations have been made via a multi-domain 
spectral method |3Jj] . We have developed a code to imple- 
ment this method by using the C++ library Lorene |3lj . 
The key advantage of the spectral method is that the re- 
quired number of grid points to obtain a sufficiently high 
accuracy is quite small compared to the grid points in 
finite differencing, and the accuracy is guaranteed up to 
the round-off error in principle. Since this method is only 
applicable to smooth functions, we treat the discontinuity 
at the surface of the star by splitting the computational 
domain. Note that we introduce three computational do- 
mains to cover the space. The innermost domain covers 
the whole star, while the outermost one is compactificd 
which allows to cover the space up to spatial infinity. We 
also use surface fitting method to split the domain at the 
surface of the star [30J . This method works perfect until 
the equatorial surface of the star has a cusp, which hap- 
pens when the uniformly rotating star approaches to the 
mass shedding limit, or the star is highly deformed from 
the sphere due to differential rotation. 



B. Iterative evolution approach 

We follow the iterative evolution approach 0, to 
investigate the viscosity driven instability in rotating rel- 
ativistic stars. We particularly focus on the effect of rel- 
ativistic gravitation for compressible fluids. The phys- 
ical viewpoint of this approach is only shown in New- 
tonian incompressible star that to study the transition 
from a uniformly rotating axisymmetric body (Maclau- 
rin spheroid) to a nonaxisymmetric body (Jacobi ellip- 
soid). According to Christodoulou et al. [32], the above 
deformation process is driven by viscosity, since it only 
varies the circulation but keeps the other two conserved 
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| Initial guess: Spherically symmetric starsj 








■* |-Guess of next matter configuration | 


1 . Impose rotation at the given iteration step 

2. Start varying the angular velocity 

to the required value at the given iteration step 

3. Stop varying the angular velocity 

at the required value at the given iteration step 





Solve gravitational field equations 



| Check the convergence of the norm of the enthalpy | 

| YES 

| Equilibrium configuration] 

| < — Impose a bar mode perturbation 
| Check the signature of q/q~| 
positive^ ^ ^ neg ative 

| Unstable | | Stable | 

FIG. 1: Sketch of the computational procedure to determine the stability of rotating stars. 



quantities, total energy and angular momentum, in the 
Newtonian incompressible star. The computational view- 
point of this approach is that instead of performing the 
time evolution of the star to investigate the stability of 
the star, we treat the iterative number as evolutional 
time and determine the stability of the star. The ad- 
vantage of this approach is that there is no restriction 
to the evolutional timestep even in a high compactness 
star. Note that it is entirely difficult to study the secular 
instability in the explicit evolution scheme in full general 
relativistic hydrodynamics, since the restriction from the 
Courant timestep scales as ~ {M / R)^ 1 / 2 . The uncer- 
tain issue in this approach is that whether one can treat 
iterative number as evolutional time, and there is no re- 
lationship between each other in a mathematical sense. 
However there is a correspondence: In the Newtonian 
incompressible star, Gondek-Rosihska and Gourgoulhon 
[l7| investigate the difference between the exact critical 
value on the bifurcation point (e.g. |f|) and find that it 
is within the round-off error. 

Our computational study of the iterative evolution ap- 
proach is divided into two stages; construction of a rotat- 
ing equilibrium star and the determination of the viscos- 
ity driven bar mode stability of a rotating equilibrium 
star. To construct a rotating equilibrium star, we first 
construct a spherical star with given parameters of T and 
ff mal , as an initial guess of the metric components and 
the matter profile. Next we solve gravitational field equa- 
tions and determine the new matter profile for the next 
iteration step. During the iteration we impose rotation 
at the given iteration step, start varying the angular ve- 
locity to the given required value at the given iteration 
step, stop varying the angular velocity at the given re- 
quired value at the given iteration step. We stop our 



iteration cycle when the relative error of the enthalpy 
norm between the previous step and the current step is 
within 1.0 x 10~ 7 . We check the relativistic virial identity 
and the identity for all the equilibrium configuration and 
find that the relative errors of GRV2 and GRV3 are <■ 
several xl0~ 4 . 

To determine the stability of rotating relativistic star 
driven by viscosity, we follow all the computational pro- 
cedure to construct a rotating equilibrium star until we 
reach the relative error of the enthalpy norm at 1.5 x 10~ 7 . 
At this iteration step, we put the following m = 2 per- 
turbation in the logarithmic lapse to enhance the growth 
of the bar mode instability as 

v = v cq (l + e amp sin 2 6 cos 2ip), (2.39) 

where v cq is the logarithmic lapse in the equilibrium, e amp 
is the amplitude of the perturbation. We diagnose the 
maximum logarithmic lapse of the m = 2 coefficients in 
terms of mode decomposition as 

g = max|j> 2 |, (2.40) 

where 

oo 

v = J2 v m e mV - (2-41) 

m—0 

We also define the logarithmic derivative of q in the iter- 
ation step Mi as 



where denotes q at the iteration step A/i. We deter- 
mine the stability of the star in terms of viscosity driven 
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FIG. 2: Diagnostic q/q as a function of iteration numbers M 
for sixteen different uniformly rotating stars of Y = 2.3 and 
M/R = 0.01. Solid, dash-dotted, dashed, dotted lines de- 
notes T/W = (0.1306, 0.1307, 0.1308, 0.1309) for {Me, e amp ) = 
(17,1.00 x 10~ 3 ), T/W = (0.1309,0.1312,0.1313,0.1314) 
for (7V e ,e amp ) = (17,1.00 x HT 4 ), T/W = 
(0.1309,0.1312,0.1313,0.1314) for (7V 9 ,£ amp ) 
(17, 1.00 x 10~ 5 ), and T/W = (0.1306, 0.1307, 0.1308, 0.1310) 
for (A/s,e amp ) = (25, 1.00 x 10~ 5 ), respectively. 



instability as follows. When the diagnostic q grows expo- 
nentially after we impose a bar mode perturbation in the 
logarithmic lapse, we conclude that the star is unstable. 
On the other hand when the diagnostic decays after we 
put a perturbation, the star is stable. More precisely, 
we monitor the derivative of q after we impose a bar- 
mode perturbation. We conclude that the equilibrium 
star is unstable when the q/q settles down to a positive 
constant value, while stable to a negative constant value. 
Note that the existence of the plateau in q/q after the 
several iteration steps once we put a perturbation con- 
firms us that we are in the linear perturbation regime, 
and therefore guarantees our choice of the perturbation 
amplitude we imposed (e am p)- Finally we determine the 
critical value of T/W as the minimum one in the unstable 
branch. We also confirm our argument in all equilibrium 
stars that there is a continuous transition between stable 
and unstable stars as a function of T/W in our model. 
We summarize our computational procedure in Fig. ^ 



III. UNIFORMLY ROTATING STARS 

Before studying the viscosity driven instability in uni- 
formly rotating stars, we examine the dependence of the 
perturbation amplitude £ amp and the collocation points 
on the critical value of viscosity driven bar mode insta- 
bility (T/W) cr t- Note that we introduce three domains 
to cover the whole space, and each domain has a rela- 
tionship of M r = 2Mb — 1 and M v = 4, where M r , Mg, M^ 
represents the collocation points for the radial direction, 
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FIG. 3: Stability of uniformly rotating stars of V = 2.3, 
M/R = 0.01 for four different parameters. Circle (open) and 
circle (filled) denotes stable and unstable to viscosity driven 
bar mode perturbation respectively. We fix the compactness 
for each parameter up to four digits. Note that there is a 
monotonic transition from stable to unstable when increasing 
T/W. 



the polar angle, the azimuthal angle, respectively. 

First we vary the amplitude e am p in the range between 
1.00 x 10~ 3 - 1.00 x 10 -5 . We show the diagnostic q/q in 
Fig- El The general picture of q/q in our computation is 
composed of three stages: (1) q/q — (2) Sudden change 
in q/q at a certain iteration step (3) Continuous smooth 
function. (1) represents the stage of constructing the 
axisymmctric rotating equilibrium configuration. Since 
we solve the equations in the axisymmetric spacetime, 
there is no non-axisymmetric contribution in this stage 
and therefore both q and q/q are up to numerical error. 
(2) represents a reaction due to a sudden imposition of 
a bar-mode perturbation in the logarithmic lapse. The 
quantity q/q should drastically change at this iteration 
step. (3) represents the post-perturbation stage, which 
determines the stability of the equilibrium star. When 
q/q takes a positive value after a certain iteration step 
we conclude that the equilibrium star is unstable, while 
takes a negative value stable. Note that it is in the linear 
perturbation regime when q/q takes a constant value af- 
ter the perturbation. Therefore the solid and dash-dotted 
lines in Fig. El are stable while dash and dotted lines are 
unstable (plotted in Fig.|3J). The diagnostic q/q suggests 
us to impose the amplitude below 1.00 x 10~ 4 , since there 
is no plateau for the case of Mg = 17, £ amp = 1.00 x 10 -3 
in the stable star (Fig.EJ). We also check whether we have 
a continuous transition from the stable star to the unsta- 
ble one in terms of T/W (Fig. [5J, and confirmed that 
there exists a minimum T/W in the unstable stars. We 
determine the minimum value as (T/W) cr t. The mono- 
tonic increase of the final value of q/q as increasing T/W 
(Fig. E| also supports the previous statement. We find 
in Table [I] that the critical value of T/W depends on the 
choice of e amp for only 0.4%, which means that we are in 
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FIG. 4: Diagnostic q/q as a function of iteration steps M for sixteen different uniformly rotating stars. Solid and dashed lines 
denotes the unstable and stable stars respectively. Note that the T/W for each stable star in the same compactness is 0.0001 
lower than the critical value of that of an unstable star (see Table ITU . 



the high convergence level of the choice of e amp . 

Next we show our result of two different choices of the 
collocation points Me — 17 and 25 in Table [I] and find 
that the critical T/W depends on the choice of M r and 
Me only for 0.4%. This means that we are also in the 
high convergence level in the choice of collocation points. 
Therefore we briefly estimate that our accuracy level of 
the critical value of T/W is < 1%. Hereafter we choose 
the parameter sets e a m P = 1.00 x 10 -5 and Me — 17 to 
determine the critical value of viscosity driven instability 
in uniformly rotating stars. 

TABLE I: Dependence of critical T/W on the amplitude e a m P 
and on the collocation points in uniformly rotating stars of 
T = 2.3, M/R = 0.01. 



£amp 




Me 


(T/W0crt 


1.00 x 10" 


3 


17 


0.1308 


1.00 x 10" 


4 


17 


0.1313 


1.00 x 10" 


5 


17 


0.1313 


1.00 x 10" 


5 


25 


0.1308 



After we determine our choice of e amp and the collo- 
cation points, we study the critical value of T/W of the 
viscosity driven instability in uniformly rotating stars. 
We show the diagnostic of the two closest star to the 
critical value of T/W for sixteen different parameters 
in Fig. Note that there is a clear plateau after we 
put a perturbation, and therefore we are in the linear 
perturbation regime. We also study the stability of the 
stars in terms of T/W (Fig.|SJ and determine the critical 
values of instability, which are summarized in Table [H] 
and in Fig. |SJ We find that the relativistic gravitation 
stabilizes the star from viscosity driven instability, and 
that the critical value of T/W for each compactness is 
almost insensitive to the polytropic T of the equation of 
state. The Newtonian compressible calculation has been 
performed in Fig. 3 of Bonazzola, Frieben, and Gour- 
goulhon [HJ that the critical value of T/W is - 0.134 
which is not so sensitive to the variation of the poly- 
tropic r. Our computational results (Fig. ||JJ) shows that 
the critical T/W is - 0.137 for M/R = 0.01 , - 0.145 
for M/R = 0.05, - 0.150 for M/R = 0.1, - 0.157 for 
M/R = 0.15, respectively. The critical value of T/W 
monotonically increases when increasing the compactness 
of the star, which means that the relativistic gravitation 
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FIG. 5: Stability of uniformly rotating stars for sixteen different parameters. Circle (open) and circle (filled) denotes stable 
and unstable to viscosity driven bar mode perturbation respectively. We fix the compactness for each parameter up to four 
digits. Note that there is a monotonic transition from stable to unstable when increasing T/W. 



stabilizes the viscosity driven instability. Note that the 
r below the smallest T in each compactness of the star 
plotted in Fig. represents that the star is stable up to 
the mass-shedding limit. 



IV. DIFFERENTIALLY ROTATING STARS 

Next we follow the same approach to study the critical 
value of viscosity driven bar-mode instability in differen- 
tially rotating stars. We study two cases based on the 
variation of rotation profile during the iteration for dif- 
ferentially rotating stars. One is that we fix the rotation 
profile throughout the iteration. The physical timestep 
which corresponds to the iteration step is much shorter 
than the dynamical time in this case since the different 
fragments of the fluid in terms of a cylindrical radius 
move at different azimuthal speeds, and therefore the 
trace shows a spiral structure. The other is to change 
the rotation profile throughout the evolution. Since we 
mimic the model that the viscosity changes the rotation 
profile, the total angular momentum is approximately 



conserved throughout the process. We use the same col- 
location points as in the case of uniformly rotating stars. 

First we show the case of fixed rotation profile through- 
out the evolution. The diagnostic in Fig. [7\ shows that 
there exists a critical value of viscosity driven instabil- 
ity. Note that the plateau at the late stage of the itera- 
tion clearly shows that we are in the linear perturbation 
regime. We also monitor the stability of differentially ro- 
tating stars in terms of different T/W in Fig. [S] We find 
a monotonic transition from a stable star to an unstable 
star as a function of T/W, which guarantees the determi- 
nation of critical value of T/W as a minimum T/W in the 
unstable stars. We summarize our finding in Table ITTT1 

We find the following two issues in the critical value 
of T/W. One is that relativistic gravitation also stabi- 
lizes differentially rotating stars from the viscosity driven 
instability. The above statement also holds in uniformly 
rotating incompressible relativistic stars and in uniformly 
rotating compressible relativistic stars, and therefore we 
find that this statement is quite general one. The other 
is that differential rotation also stabilize the star from 
the viscosity driven instability. The critical value of 
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TABLE II: Critical value of T/W of viscosity driven instability in uniformly rotating relativistic stars 



r " 


Rp /Re b 


TT C 

max 




R d 


M e 


J* 


{T/W)*t 


M/R 


GRV2 




GRV3 




2.30 


0.5276 


9.199 x 10" 


-3 


0.9323 


0.009322 


0.0002082 


0.1313 


0.009999 


6.44 x 10" 


5 


-1.54 x 10 


-4 


2.40 


0.5446 


8.465 x 10" 


-3 


0.7518 


0.007518 


0.0001400 


0.1336 


0.01000 


3.22 x 10" 


5 


-9.60 x 10 


-5 


2.50 


0.5526 


8.002 x 10" 


-3 


0.6322 


0.006322 


0.0001008 


0.1350 


0.01000 


-3.78 x 10 


-5 


1.04 x 10" 


4 


2.60 


0.5482 


7.662 x 10" 


-3 


0.5528 


0.005528 


7.904 x 10~ 5 


0.1391 


0.01000 


5.04 x 10" 


5 


-2.14 x 10 


-5 


2.70 


0.5590 


7.368 x 10" 


-3 


0.4825 


0.004825 


6.057 x 10~ 5 


0.1385 


0.01000 


7.99 x 10" 


5 


-1.73 x 10 


-4 


2.80 


0.5510 


7.130 x 10" 


-3 


0.4352 


0.04352 


5.023 x 10"" 5 


0.1420 


0.01000 


-1.24 x 10 


-5 


8.59 x 10" 


5 


2.90 


0.5663 


6.933 x 10" 


-3 


0.3879 


0.003879 


3.966 x 10" 5 


0.1389 


0.01000 


5.85 x 10" 


5 


-1.52 x 10 


-4 


2.50 


0.5312 


4.542 x 10" 


2 


0.8140 


0.04071 


0.001986 


0.1409 


0.05001 


4.72 x 10" 


5 


-1.30 x 10 


-4 


2.60 


0.5368 


4.291 x 10" 


-2 


0.7267 


0.03633 


0.0001616 


0.1434 


0.05000 


5.04 x 10" 


5 


-4.94 x 10 


-5 


2.70 


0.5486 


4.090 x 10" 


-2 


0.6513 


0.03256 


0.001306 


0.1426 


0.05000 


7.01 x 10" 


5 


-1.64 x 10 


-4 


2.80 


0.5380 


3.951 x 10" 


-2 


0.6089 


0.03044 


0.001171 


0.1483 


0.05000 


-1.02 x 10 


-4 


1.34 x 10" 


4 


2.90 


0.5471 


3.820 x 10" 


-2 


0.5586 


0.02793 


0.0009868 


0.1469 


0.05000 


4.99 x 10" 


5 


-1.23 x 10 


-4 


2.70 


0.5241 


9.710 x 10" 


2 


0.7320 


0.07319 


0.005049 


0.1501 


0.09999 


7.93 x 10" 


r, 


-7.16 x 10 


-5 


2.80 


0.5351 


9.197 x 10" 


-2 


0.6734 


0.06735 


0.004319 


0.1506 


0.1000 


7.23 x 10" 


5 


-1.86 x 10 


-4 


2.90 


0.5423 


8.824 x 10" 


-2 


0.6274 


0.06274 


0.003773 


0.1507 


0.1000 


8.09 x 10" 


5 


-3.63 x 10 


—5 


2.90 


0.5242 


1.608 x 10" 


1 


0.6634 


0.09957 


0.008380 


0.1567 


0.1501 


3.16 x 10" 


5 


-6.99 x 10 


-5 



T: Adiabatic index of the equation of state 

b R p /R e '. Ratio of the polar proper radius to the equatorial proper 
radius 

c Hmax'- Maximum enthalpy 
d R: Circumferential radius 
e M: Gravitational mass 
fj: Total angular momentum 
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FIG. 6: Critical value of T/W as a function of Y for four dif- 
ferent compactness of uniformly rotating stars (see Table UTt . 
Circle (open), circle (filled), square (open), and square (filled) 
denotes the compactness (M/R) of 0.01, 0.05, 0.1, and 0.15, 
respectively. The star whose adiabatic index is ri ow — 0.1, 
where Fi ow is the lowest V of the unstable star in each com- 
pactness in the figure, is stable. 



T/W is 0.13 ~ 0.16 for a uniformly rotating star, while 
0.18 ~ 0.25 for a differentially rotating stars with mod- 
erate degree of differential rotation. 

Next we study the variation of the rotation profile since 
viscosity also takes a significant role to change it. The 
timescale between the change of angular velocity profile 



and the growth of the viscosity driven bar-mode are in 
the same order in low viscosity approximation, we should 
check whether our result is significantly affected by the 
change of the rotation profile of the star. In order to 
mimic this idea, we vary the parameter of the degree of 
differential rotation slightly after we impose a perturba- 
tion in the following manner: 



/I" 1 — A 

^rot — ^rot 



1 (eq) 



[l-e rot (AA-A^ p tb)], 



(4.1) 



where e ro t is the degree of the variation of the rotation 
profile we set to 1.0 x 10~ 4 , M the iteration number, 
A/'ptb the iteration number we impose the perturbation. 
Note that viscosity changes the rotation profile to the 
uniform one, we put negative sign in front of £ ro t. Since 
the viscosity only affects the local interaction between 
the each fluid components, the total angular momentum 
is conserved even the viscosity takes the role. Therefore 
we also vary the central angular velocity il c in a following 
manner 



^c = ^ Cq) [l 



; (A/--A/- P tb)], 



(4.2) 



where £ omg is the degree of the variation of the central 
angular velocity in order to maintain the total angular 
momentum approximately constant. Note that we put 
negative sign in front of e omg to play an appropriate role 
of the viscosity. In practice we vary e omg in 0.1 steps 
and choose the one that changes the angular momentum 
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FIG. 7: Diagnostic q/q as a function of iteration steps Af for five different differentially rotating stars. Solid and dashed lines 
denotes the unstable and stable stars respectively. Note that the T/W for each stable star in the same compactness is O.OOOf 
lower than the critical value of that of an unstable star (see Table ITTTt . 



TABLE III: Critical value of viscosity driven instability in differentially rotating relativistic stars. We choose V — 2 polytropic 
equation of state and the degree of differential rotation as A — 1. 



Rp 1 Re 






R 


M 


J 


(T/W)ctt 


M/R 


GRV2 


GRV3 


0.4458 


7.594 x 10" 


3 


1.978 


0.01978 


0.01165 


0.1828 


0.01000 


-2.03 x 10" 6 


2.28 x 10~ 6 


0.3985 


3.830 x 10" 


2 


1.894 


0.09468 


0.01315 


0.1999 


0.05000 


3.92 x 10" 6 


4.92 x 10" 4 


0.3820 


7.492 x 10" 


2 


1.758 


0.1758 


0.03565 


0.2186 


0.1000 


3.88 x 10" 7 


-4.52 x 10~ 6 


0.3457 


1.116 x 10" 


1 


1.609 


0.2414 


0.06025 


0.2354 


0.1500 


-1.93 x 10" 5 


-9.03 x 10" 6 


0.2982 


1.581 x 10" 


1 


1.455 


0.2910 


0.08210 


0.2496 


0.2000 


-3.63 x 10 -6 


-1.55 x 10" 5 
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FIG. 8: Stability analysis of differentially rotating stars for 
five different parameters. Circle (open) and circle (filled) de- 
notes stable and unstable to viscosity driven bar mode per- 
turbation respectively. We fix the compactness for each pa- 
rameter up to four digits. Note that there is a monotonic 
transition from stable to unstable when increasing T/W. 



minimum at the given e rot . For each case the relative 
change of the total angular momentum after we impose 
the perturbation is in the order of < 10~ 5 . We show 
the relationship between the relative change of the total 
angular momentum and e om g for several cases in Fig. 

Taking into account of the change of rotational profile, 
we show our numerical results for the critical values in 
Fig. ^1] We find that all the stars around the critical 
T/W determined for a fixed rotation profile becomes un- 
stable. In fact, the stage of q/q that corresponds to the 
plateau in Fig. [7| increases in Fig. ^| Therefore we esti- 
mate the relevant two timescales, growth of the bar mode 
due to viscosity and the variation time of the rotational 
profile due to viscosity, and compare them to discuss the 
condition to induce viscosity driven instability. 

If we assume that the bar grows exponentially through- 
out the iteration, the growth timescale of the bar is 
{q/q) -1 - T ne value of q/q at the plateau in Fig.[7]repre- 
sent the growth timescale. On the other hand, the change 
of the rotation profile due to viscosity changes the growth 
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FIG. 9: Relative deviation of the total angular momentum between two successive steps in the iteration process, as a function 
of the step number TV. Solid, dashed, dotted, dash-dotted lines represent £ omg = (0.5, 0.6, 0.7, 0.8) x 10~ 4 for M/R — 0.01 
and T/W = 0.1828, £ omg = (0.8,0.9,1.0,1.1) x 10~ 4 for M/R = 0.05 and T/W = 0.1999, e omg = (1.1,1.2,1.3,1.4) x 10~ 4 



for M/R = 0.1 and T/W 



(1.0,1.1,1.2,1.3) x 10" 4 for M/R = 0.15 and T/W = 0.2354, and e omg = 



(0.9, 1.0, 1.1, 1.2) x 10~ 4 for M/R = 0.2 and T/W = 0.2496, respectively. 




FIG. 10: Diagnostic q/q as a function of iteration steps Af for five different differentially rotating stars. Solid, dashed, 
dotted, and dashed line denotes T/W = (0.1825,0.1828,0.1834,0.1844) and e omg = 0.8 x 10~ 4 for M/R = 0.01, T/W = 
(0.1993,0.1996,0.1998,0.1999) and e omg = 0.9 x 10~ 4 for M/R = 0.05, T/W = (0.2180,0.2184,0.2185,0.2186) and e omg = 
1.3 x 10~ 4 for M/R = 0.1, T/W = (0.2348,0.2352,0.2353,0.2354) and e omg = 1.1 x 10~ 4 for M/R = 0.15, and T/W = 
(0.2494,0.2495,0.2496) and e omg = 1.0 x 10~ 4 for M/R = 0.2, respectively. Note that q/q is always increasing around the 
critical value of T/W in differentially rotating stars, which means the change of rotational profile due to viscosity unstabilizes 
the star. 
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timescale of the bar, and therefore the derivative of the 
q/q represents the timescale of the change of rotation 
profile. We define those two timescales as 



(<f/9) 



-1/2 



T ba 



(4.3) 



We summarize those two timescales from our computa- 
tion in Table |IV| 



TABLE IV: Two timescales in differentially rotating stars 
extracted from our numerical results 



M/R 


T/W 


Tang [AN] 


r bar [AN] 


0.01 


0.1825 


1.8 x 10 2 


-6.1 x 10 2 


0.01 


0.1828 


1.7 x 10 2 


5.1 x 10 4 


0.01 


0.1834 


1.3 x 10 2 


3.0 x 10 2 


0.01 


0.1844 


8.1 x 10 1 


1.1 x 10 2 


0.05 


0.1993 


1.7 x 10 2 


-1.7 x 10 3 


0.05 


0.1996 


1.6 x 10 2 


-2.7 x 10 3 


0.05 


0.1998 


1.6 x 10 2 


-8.1 x 10 3 


0.05 


0.1999 


1.6 x 10 2 


3.1 x 10 5 


0.10 


0.2180 


1.6 x 10 2 


-6.6 x 10 2 


0.10 


0.2184 


1.6 x 10 2 


-2.6 x 10 3 


0.10 


0.2185 


1.5 x 10 2 


-8.1 x 10 3 


0.10 


0.2186 


1.5 x 10 2 


8.0 x 10 3 


0.15 


0.2348 


2.0 x 10 2 


-4.8 x 10 2 


0.15 


0.2352 


2.0 x 10 2 


-1.7 x 10 3 


0.15 


0.2353 


1.9 x 10 2 


-3.4 x 10 4 


0.15 


0.2354 


1.9 x 10 2 


5.9 x 10 3 


0.20 


0.2494 


8.1 x 10 2 


-1.8 x 10 3 


0.20 


0.2495 


8.5 x 10 2 


-2.4 x 10 4 


0.20 


0.2496 


7.7 x 10 2 


2.3 x 10 3 



These two timescales can also be derived analytically 
in Newtonian gravity. The azimuthal component of the 
Newtonian Navier-Stokes equation is 



~dt 



v d 
vj 3 dvj 



dm 



(4.4) 



where v is a shear viscosity. We assume the time de- 
pendence of the angular velocity as il(t) — S! oq e -I °'* at 
the derivation from the equilibrium due to viscosity, the 
timescale to change the angular velocity profile to the 
uniform rotation is 



(zu 2 + d 2 ) 2 d 2 r 2 n c 
> — = 



8d 2 v 



8v fL — fh 



(4.5) 



where Q s is the equatorial surface angular velocity. Note 
that we adopt the j-constant rotational profile to the 
angular velocity (Eq. \2.'2'2\ in Newtonian gravity). The 
e-folding time of a compressible uniformly rotating star 
(Eq. [A20] of Ref. ) based on the Navier-Stokes equa- 
tion is 



T~ha 



Pscc 

0-Ps. 



(4.6) 



where structure constant depends on the poly- 

tropic index (Table 1 of Ref. [H; n n = 0.65345 for 
r = 2), (3 SCC a critical value of f3 for the secular bar mode 
instability, R is the equatorial radius of the star. 

Based on the analytical estimation of the two 
timescales, the timescales used in our numerical results 
should be described as 



T ba 



-■ore ) 



-org 



/9-/3L 



(4.7) 
(4.8) 



We confirm that differential rotation stabilizes the star 
from viscosity driven instability. Also Eq. i|4.8[l shows 
that the timescale of the bar becomes short when the ro- 
tational profile varies due to viscosity. The critical value 
of T/W changes from the one computed in a fixed rota- 
tion profile, but the deviational ratio of T/W is roughly 
the same order as the one of rotational profile, which 
means « e re- 



V. CONCLUSIONS 

We have studied the viscosity driven instability in both 
uniformly and differentially rotating polytropic star by 
means of iterative evolution approach in general relativ- 
ity. We have focussed on the determination of the critical 
value of viscosity driven instability. 

We find that relativistic gravitation stabilizes the star 
from the viscosity driven instability, with respect to New- 
tonian gravitation. Also the critical value is not sensitive 
to the stiffness of polytropic equation of state for a given 
compactness of the star. In a previous study devoted to 
compressible stars, Bonazzola, Frieben, and Gourgoulhon 
[l6l | showed that relativistic gravitation does stabilize the 
uniformly rotating polytropic stars by investigating the 
mass-shedding sequence. Our study shows the concrete 
value of (T/W) CI t in the case of a uniformly rotating star. 
Also we have improved the numerical code with respect 
to the study 0] by making use of surface-fitted coordi- 
nates. 

Besides we have found that differential rotation also 
stabilizes the star from the viscosity driven instability. 
If we fix the compactness of the star, we find a signifi- 
cant increase of the critical value of T/W, which supports 
the above statement. We also confirm this statement by 
changing the rotation profile due to viscosity and find 
that the differential rotation still stabilizes the threshold 
of viscosity driven instability. However to confirm the 
statement in differentially rotating stars, the study of 
other approaches such as implicit hydrodynamical evolu- 
tion or eigenmode analysis should be necessary and help- 
ful. 

Finally let us mention the characteristic amplitude and 
frequency of gravitational waves emitted throughout the 
viscosity driven secular bar mode instability, which pro- 
duces quasi-periodic gravitational waves detectable in 
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ground based interferometers. The characteristic ampli- 
tude h of gravitational waves estimated from the evo- 



where d is the distance to the source and the charac- 
teristic frequency / = fl/n is / > 1000 [Hz], depending 
on T/W of the star. Note that the frequency increases 
throughout this process. Although the frequency regime 
of the source is slightly higher than the best sensitivity 
regime of the ground based detectors to follow all the de- 
formation process, we may have a chance to detect them 
when it happens in the Virgo cluster, for example. 
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